
	cd(datadir)
	
	# Packages
		# To install: using Pkg
		# Then Pkg.add("package_name")
		# NOTE: I'm not sure we need all of these packages
		
		using DataFrames
		using StatsBase
		using Optim
		using Calculus
		using Distributions
		using NLsolve
		using ForwardDiff
		using Random
		using DelimitedFiles
		using CSV
		using LinearAlgebra
		using LineSearches
	
	Random.seed!(6);
	
##### Main functions that calculate terms in FOCs
	
	cd(codedir)
	include("estimation.functions.jl");	
	include("counterfactual.functions.jl");	
	cd(datadir)
	
##### Set Up Counterfactual Scenario
	
	scenario_definitions = CSV.read("scenario_definitions_inertia.csv", DataFrame, header=true);
	scenarios = scenario_definitions[!,:Scenario];
	
	s_index = findall(scenarios .== scenario)[1];
	nora_flag = scenario_definitions[s_index,:nora_flag];
	mand_repeal_flag = scenario_definitions[s_index,:mand_repeal_flag];
	voucher_flag = scenario_definitions[s_index,:voucher_flag];
	no_inertia_flag = scenario_definitions[s_index,:no_inertia_flag]; 
	pc_flag = scenario_definitions[s_index,:pc_flag];
	no_uninsured_flag = scenario_definitions[s_index,:no_uninsured_flag];
	except_network_flag = scenario_definitions[s_index,:except_network_flag];
	nochurn_flag = scenario_definitions[s_index,:nochurn_flag];
	network_flag = scenario_definitions[s_index,:network_flag];
	eq_robust_flag = scenario_definitions[s_index,:eq_robust_flag];
	inattention_flag = scenario_definitions[s_index,:inattention_flag];
	search_flag = scenario_definitions[s_index,:search_flag];
	inattention = inattention_flag;
	
	if create_standard_errors_flag
		premium_file = string(year,scenario,"0_premium_output_intvars_avint_",sedraw,".csv");
		output_file = string(year,scenario,"0_supply_output_intvars_avint_",sedraw,".csv");
	elseif (choice_seq_flag & int_rs_variables & add_av_interactions)
		premium_file = string(year,scenario,"0_premium_output_intvars_seq_avint.csv");
		output_file = string(year,scenario,"0_supply_output_intvars_seq_avint.csv");
	elseif (choice_seq_flag & int_rs_variables)
		premium_file = string(year,scenario,"0_premium_output_intvars_seq.csv");
		output_file = string(year,scenario,"0_supply_output_intvars_seq.csv");
	elseif (all_rs_variables & add_av_interactions)
		premium_file = string(year,scenario,"0_premium_output_allvars_avint.csv");
		output_file = string(year,scenario,"0_supply_output_allvars_avint.csv");
	elseif (int_rs_variables & add_av_interactions)
		premium_file = string(year,scenario,"0_premium_output_intvars_avint.csv");
		output_file = string(year,scenario,"0_supply_output_intvars_avint.csv");
	elseif int_rs_variables
		premium_file = string(year,scenario,"0_premium_output_intvars.csv");
		output_file = string(year,scenario,"0_supply_output_intvars.csv");
	elseif hisp_flag
		premium_file = string(year,scenario,"0_premium_output_hisp.csv");
		output_file = string(year,scenario,"0_supply_output_hisp.csv");
	elseif agegender_flag
		premium_file = string(year,scenario,"0_premium_output_agegender.csv");
		output_file = string(year,scenario,"0_supply_output_agegender.csv");
	elseif old_flag
		premium_file = string(year,scenario,"0_premium_output_oldflag.csv");
		output_file = string(year,scenario,"0_supply_output_oldflag.csv");
	elseif zero_resid_flag & choice_seq_flag
		premium_file = string(year,scenario,"0_premium_output.csv");
		output_file = string(year,scenario,"0_supply_output.csv");
	elseif zero_resid_flag
		premium_file = string(year,scenario,"0_premium_output_noseq.csv");
		output_file = string(year,scenario,"0_supply_output_noseq.csv");
	elseif choice_seq_flag
		premium_file = string(year,scenario,"_premium_output.csv");
		output_file = string(year,scenario,"_supply_output.csv");
	else
		premium_file = string(year,scenario,"_premium_output_noseq.csv");
		output_file = string(year,scenario,"_supply_output_noseq.csv");
	end
	
##### Load Data		
	
	if network_flag
		data = CSV.read("ca_julia_supply_data_JUN182021_net.csv", DataFrame, header=true); 
		households = CSV.read("ca_household_characteristics_JUN182021_net.csv", DataFrame, header = true);
		unique_choice_sets = CSV.read("unique_choice_sets_JUN182021_net.csv", DataFrame,header = true);
	elseif nochurn_flag
		data = CSV.read("ca_julia_supply_data_JUN262020_churn.csv", DataFrame, header=true); 
		households = CSV.read("ca_household_characteristics_AUG032019_churn.csv", DataFrame, header = true);
		unique_choice_sets = CSV.read("unique_choice_sets_MAY112021_churn.csv", DataFrame,header = true);
		#foc_residuals = convert(Array,CSV.read(string(year,"_foc_residuals.csv"), DataFrame,header = false)[!,1]);
	else
		data = CSV.read("ca_julia_supply_data_JUN262020_learn.csv", DataFrame, header=true); 
		households = CSV.read("ca_household_characteristics_AUG032019_small.csv", DataFrame, header = true);
		unique_choice_sets = CSV.read("unique_choice_sets_MAY112021.csv", DataFrame,header = true);	
	end
	
	if network_flag
		
		# Plan characteristic data w/o CSR plans
		plans = CSV.read("ca_plan_year_JUN182021_net.csv", DataFrame, header = true);
		plans_pmt = CSV.read("ca_plan_market_year_JUN182021_net.csv", DataFrame,header = true);
		
		# Demand parameter estimates
		param_estimates = CSV.read("demand_parameter_estimates.csv", DataFrame, header=true); # parameters from demand estimation
		param_estimates = dropmissing(param_estimates,:Network);
		min_beta = param_estimates[!,:Network];
		
		# Claims and risk score parameter estimates
		claims_parameter_estimates = CSV.read("claims_estimates_V11_net.csv",DataFrame,header=true);
		rs_parameter_estimates = CSV.read(string("rs_coefficients_V11_inertia_net.csv"),DataFrame,header=true);
		
	elseif eq_robust_flag
	
		# Plan characteristic data w/o CSR plans
		plans = CSV.read("ca_plan_year_JUN262020_learn.csv", DataFrame, header = true);
		plans_pmt = CSV.read("ca_plan_market_year_JUN262020_learn.csv",DataFrame,header = true);
		
		# Demand parameter estimates
		param_estimates = CSV.read("demand_parameter_estimates.csv",DataFrame,header=true); # parameters from demand estimation
		param_estimates = dropmissing(param_estimates,:Robust);
		min_beta = param_estimates[!,:Robust];
		
		# Claims and risk score parameter estimates
		claims_parameter_estimates = CSV.read("claims_estimates_V11_robust.csv",DataFrame,header=true);
		rs_parameter_estimates = CSV.read("rs_coefficients_robust.csv",DataFrame,header=true);
		
	else
		
		# Plan characteristic data w/o CSR plans
		plans = CSV.read("ca_plan_year_JUN262020_learn.csv", DataFrame, header = true);
		plans_pmt = CSV.read("ca_plan_market_year_JUN262020_learn.csv", DataFrame,header = true);
		
		# Demand parameter estimates
		if create_standard_errors_flag
			Random.seed!(sedraw);
			param_estimates = CSV.read("base_demand_parameter_estimates.csv", DataFrame, header=true); # parameters from demand estimation
			min_beta = param_estimates[!,:Coef] .+ randn(size(param_estimates)[1]) .* param_estimates[!,:SE] .* param_estimates[!,:Sensitivity];
			Random.seed!(6);
		else
			param_estimates = CSV.read("demand_parameter_estimates.csv", DataFrame, header=true); # parameters from demand estimation
			if (choice_seq_flag & add_av_interactions)
				param_estimates = dropmissing(param_estimates,:AV_Seq_CF_Random_IMFE);
				min_beta = param_estimates[!,:AV_Seq_CF_Random_IMFE];
			elseif choice_seq_flag
				param_estimates = dropmissing(param_estimates,:Seq_CF_Random_IMFE);
				min_beta = param_estimates[!,:Seq_CF_Random_IMFE];
			elseif add_av_interactions
				param_estimates = dropmissing(param_estimates,:AV_CF_Random_IMFE);
				min_beta = param_estimates[!,:AV_CF_Random_IMFE];
			else
				param_estimates = dropmissing(param_estimates,:CF_Random_IMFE);
				min_beta = param_estimates[!,:CF_Random_IMFE];
			end
		end
		
		# Claims and risk score parameter estimates
		if (old_flag & av_rs_flag)
			claims_parameter_estimates = CSV.read("claims_estimates_old_av.csv",DataFrame,header=true);
		elseif old_flag
			claims_parameter_estimates = CSV.read("claims_estimates_old.csv",DataFrame,header=true);
		elseif (int_rs_variables & av_rs_flag)
			claims_parameter_estimates = CSV.read("claims_estimates_int_av.csv",DataFrame,header=true);
		elseif (choice_seq_flag & int_rs_variables & add_av_interactions)
			claims_parameter_estimates = CSV.read("claims_estimates_int_seq_avint.csv",DataFrame,header=true);
		elseif (choice_seq_flag & int_rs_variables)
			claims_parameter_estimates = CSV.read("claims_estimates_int_seq.csv",DataFrame,header=true);
		elseif (all_rs_variables & add_av_interactions)
			claims_parameter_estimates = CSV.read("claims_estimates_full_avint.csv",DataFrame,header=true);
		elseif (int_rs_variables & add_av_interactions)
			claims_parameter_estimates = CSV.read("claims_estimates_int_avint.csv",DataFrame,header=true);
		elseif int_rs_variables
			claims_parameter_estimates = CSV.read("claims_estimates_int.csv",DataFrame,header=true);
		elseif hisp_flag
			claims_parameter_estimates = CSV.read("claims_estimates_hisp.csv",DataFrame,header=true);
		elseif agegender_flag
			claims_parameter_estimates = CSV.read("claims_estimates_agegender.csv",DataFrame,header=true);
		elseif (all_rs_variables & choice_seq_flag & av_rs_flag)
			claims_parameter_estimates = CSV.read("claims_estimates_full_seq.csv",DataFrame,header=true);
		elseif (all_rs_variables & av_rs_flag)
			claims_parameter_estimates = CSV.read("claims_estimates_full_av.csv",DataFrame,header=true);
		elseif all_rs_variables
			claims_parameter_estimates = CSV.read("claims_estimates_full.csv",DataFrame,header=true);
		end
		
		if (old_flag & av_rs_flag)
			rs_parameter_estimates = CSV.read(string("rs_coefficients_old_av.csv"),DataFrame, header=true);
		elseif old_flag	
			rs_parameter_estimates = CSV.read(string("rs_coefficients_old.csv"),DataFrame, header=true);
		elseif (int_rs_variables & av_rs_flag)
			rs_parameter_estimates = CSV.read(string("rs_coefficients_intermediate_av.csv"),DataFrame, header=true);
		elseif (choice_seq_flag & int_rs_variables & add_av_interactions)
			rs_parameter_estimates = CSV.read(string("rs_coefficients_intermediate_seq_avint.csv"),DataFrame, header=true);
		elseif (choice_seq_flag & int_rs_variables)
			rs_parameter_estimates = CSV.read(string("rs_coefficients_intermediate_seq.csv"),DataFrame, header=true);
		elseif (all_rs_variables & add_av_interactions)
			rs_parameter_estimates = CSV.read(string("rs_coefficients_full_avint.csv"),DataFrame, header=true);
		elseif (int_rs_variables & add_av_interactions)
			rs_parameter_estimates = CSV.read(string("rs_coefficients_intermediate_avint.csv"),DataFrame, header=true);
		elseif int_rs_variables
			rs_parameter_estimates = CSV.read(string("rs_coefficients_intermediate.csv"),DataFrame, header=true);
		elseif hisp_flag
			rs_parameter_estimates = CSV.read(string("rs_coefficients_intermediate_hisp.csv"),DataFrame, header=true);
		elseif agegender_flag
			rs_parameter_estimates = CSV.read(string("rs_coefficients_agegender.csv"),DataFrame, header=true);
		elseif (all_rs_variables & choice_seq_flag & av_rs_flag)
			rs_parameter_estimates = CSV.read(string("rs_coefficients_av.csv"),DataFrame, header=true);
		elseif (all_rs_variables & av_rs_flag)
			rs_parameter_estimates = CSV.read(string("rs_coefficients_av_noseq.csv"),DataFrame, header=true);
		elseif all_rs_variables
			rs_parameter_estimates = CSV.read(string("rs_coefficients_full.csv"),DataFrame, header=true);
		end
		
	end
	
##### Load optimal claim and risk score parameters	
	
	all_covariates = convert(Array,param_estimates[!,:Variable]);
	if nested
		pop!(all_covariates);
	end
	
	if add_silver
		product_covariates = intersect(all_covariates,
			["Premium","AV","Silver","HMO","Anthem","Blue_Shield","Kaiser","Health_Net"]);
	else
		product_covariates = intersect(all_covariates,
			["Premium","AV","HMO","Anthem","Blue_Shield","Kaiser","Health_Net"]);
	end
	
	#theta = cat(dims=1,claims_parameter_estimates[!,2],rs_parameter_estimates[!,2]); 
		
	rating_area_variables = ["share_ra2","share_ra3","share_ra4","share_ra5","share_ra6","share_ra7",  
				"share_ra8","share_ra9","share_ra10","share_ra11","share_ra12","share_ra13", 
				"share_ra14","share_ra15","share_ra16","share_ra17","share_ra18","share_ra19"];
	claims_moments_covariates = cat(dims=1,["intercept","HMO","log_risk_score","trend",
		"Anthem","Blue_Shield","Kaiser","Health_Net"],rating_area_variables);
	all_mkts = collect(2:1:19);
	
	if all_rs_variables
		risk_score_share_variables = cat(dims=1,["share_0to34","share_35to54","share_250to400","share_gt400",
			"share_male","share_family","share_Asian","share_Black","share_Hispanic","share_Other_Race"]);
	elseif int_rs_variables
		risk_score_share_variables = cat(dims=1,["share_0to34","share_male","share_family","share_minority"]);
	elseif hisp_flag
		risk_score_share_variables = cat(dims=1,["share_0to34","share_male","share_family","share_Hispanic"]);
	elseif agegender_flag
		risk_score_share_variables = cat(dims=1,["share_0to34","share_male"]);
	elseif old_flag
		risk_score_share_variables = cat(dims=1,["share_18to54","share_Hispanic"]);
	end
	
	if av_rs_flag
		risk_score_covariates = cat(dims=1,"intercept","AV",risk_score_share_variables);
	else
		risk_score_covariates = cat(dims=1,"intercept","Silver","Gold","Platinum",risk_score_share_variables);
	end
	
	premium_indices = findall(convert(Array,param_estimates[!,:Premium_index]) .== 1);	
	previous_choice_indicator = convert(Array,param_estimates[!,:Previous_Choice_Index]);
	previous_choice_indices = findall(previous_choice_indicator .== 1);
	if nested
		pop!(previous_choice_indicator);	
	end
	if network_flag	
		previous_choice_indicator_net = convert(Array,param_estimates[!,:Previous_Choice_Index_Net]);
		previous_choice_indices_net = findall(previous_choice_indicator_net .== 1);
		if nested
			pop!(previous_choice_indicator_net);		
		end
	end
	
##### Extract sample

	# First create two choice sample
	keep_households = findall(convert(Array,households[!,:year]) .<= max_year);
	if eq_robust_flag
		keep_households = intersect(keep_households,findall(convert(Array,households[!,:year]) .>= 2016));
	end

	keep_households = intersect(keep_households,findall(households[!,:year] .== year));
	keep_choice_sets = findall(unique_choice_sets[!,:year] .== year); 
	keep_choice_sets = sort(unique(data[findall(in(keep_households),data[!,:household_number]),:choice_set_id]));
	
	# If doing forward simulation, save previous year
	if year > 2014
		keep_households_ids = households[keep_households,:household_id];
		all_household_ids_in_year = findall(in(keep_households_ids),households[!,:household_id]);
		keep_households_prev = sort(intersect(all_household_ids_in_year,findall(households[!,:year] .== (year-1))));
		
		data_prev = data[findall(in(keep_households_prev),data[!,:household_number]),:];
		households_prev = households[keep_households_prev,:];
		unique_choice_sets_prev = copy(unique_choice_sets);
		plans_prev  = copy(plans);
		if (choice_seq_flag & int_rs_variables & add_av_interactions)
			premiums_prev = CSV.read(string((year-1),"_premium_output_zeroresid_intvars_seq_avint.csv"),DataFrame, header=true);
		elseif (choice_seq_flag & int_rs_variables)
			premiums_prev = CSV.read(string((year-1),"_premium_output_zeroresid_intvars_seq.csv"),DataFrame, header=true);
		elseif (all_rs_variables & add_av_interactions)
			premiums_prev = CSV.read(string((year-1),"_premium_output_zeroresid_allvars_avint.csv"),DataFrame, header=true);
		elseif (int_rs_variables & add_av_interactions)
			premiums_prev = CSV.read(string((year-1),"_premium_output_zeroresid_intvars_avint.csv"),DataFrame, header=true);
		elseif int_rs_variables	
			premiums_prev = CSV.read(string((year-1),"_premium_output_zeroresid_intvars.csv"),DataFrame, header=true);
		elseif hisp_flag
			premiums_prev = CSV.read(string((year-1),"_premium_output_zeroresid_hisp.csv"),DataFrame, header=true);
		elseif agegender_flag
			premiums_prev = CSV.read(string((year-1),"_premium_output_zeroresid_agegender.csv"),DataFrame, header=true);
		elseif old_flag
			premiums_prev = CSV.read(string((year-1),"_premium_output_zeroresid_oldflag.csv"),DataFrame, header=true);
		elseif zero_resid_flag & choice_seq_flag	
			premiums_prev = CSV.read(string((year-1),"_premium_output_zeroresid.csv"),DataFrame, header=true);
		elseif zero_resid_flag
			premiums_prev = CSV.read(string((year-1),"_premium_output_zeroresid_noseq.csv"),DataFrame, header=true);
		elseif choice_seq_flag
			premiums_prev = CSV.read(string((year-1),"_premium_output.csv"),DataFrame, header=true);
		else
			premiums_prev = CSV.read(string((year-1),"_premium_output_noseq.csv"),DataFrame, header=true);
		end
		data_prev[!,:household_number] = Int.(indexin(convert(Array,data_prev[!,:household_id]),
			convert(Array,households_prev[!,:household_id])));
		
		keep_choice_sets_prev = findall(unique_choice_sets_prev[!,:year] .== (year-1)); 
		keep_choice_sets_prev = sort(unique(data_prev[findall(in(keep_households_prev),data_prev[!,:household_number]),:choice_set_id]));
		new_choice_set_numbers_prev = zeros(Int64,size(unique_choice_sets_prev)[1]);
		new_choice_set_numbers_prev[keep_choice_sets_prev] = collect(1:1:length(keep_choice_sets_prev));
		data_prev[!,:choice_set_id] = new_choice_set_numbers_prev[data_prev[!,:choice_set_id]];
		unique_choice_sets_prev = unique_choice_sets_prev[keep_choice_sets_prev,:];
	end
	
	
	new_household_numbers = zeros(Int64,size(households)[1]);
	new_household_numbers[keep_households] = collect(1:1:length(keep_households));
	data[!,:household_number] = new_household_numbers[convert(Array,data[!,:household_number])];
	
	new_choice_set_numbers = zeros(Int64,size(unique_choice_sets)[1]);
	new_choice_set_numbers[keep_choice_sets] = collect(1:1:length(keep_choice_sets));
	data[!,:choice_set_id] = new_choice_set_numbers[data[!,:choice_set_id]];
	
	
	#data = data[findall(convert(Array,data[!,:year]) .== year),:];
	data = data[findall(data[!,:household_number] .> 0),:];
	households = households[keep_households,:];
	unique_choice_sets = unique_choice_sets[keep_choice_sets,:];
		
	uninsured_plans = convert(Array,data[!,:uninsured_plan]);
	household_numbers_data = convert(Array,data[!,:household_number]);	
	household_numbers = household_numbers_data[uninsured_plans .== 0];	
	choices = convert(Array,data[!,:choice]);
	household_weights = data[uninsured_plans .== 0,:weight];
	
	##### Remove any plans not available in sample
	
	data_pmt_names = string.(data[uninsured_plans .== 0,:plan_name],
								data[uninsured_plans .== 0,:year],
								households[household_numbers,:rating_area]);
	data_pt_names = string.(data[uninsured_plans .== 0,:plan_name],
								data[uninsured_plans .== 0,:year]);
	
	plan_names = sort(unique(data[uninsured_plans .== 0,:plan_name]));
	plan_year_names = sort(unique(data_pt_names));
	plan_pmt_names = sort(unique(data_pmt_names));
	
	plan_indicator = trues(size(plans)[1]);
	for i in eachindex(plan_indicator)
		plan_indicator[i] = in(plans[i,:pt_name],plan_year_names);
	end
	plans = plans[findall(plan_indicator),:];
	sort!(plans);

	if year > 2014
	
		data_prev[!,:household_number_link] = Int.(indexin(convert(Array,data_prev[!,:household_id]),
			convert(Array,households[!,:household_id])));
		
		uninsured_plans_prev = convert(Array,data_prev[!,:uninsured_plan]);
		household_numbers_data_prev = convert(Array,data_prev[!,:household_number]);	
		household_numbers_prev = household_numbers_data_prev[uninsured_plans_prev .== 0];	
		
		data_pt_names_prev = string.(data_prev[uninsured_plans_prev .== 0,:plan_name],
							data_prev[uninsured_plans_prev .== 0,:year]);
		plan_year_names_prev = sort(unique(data_pt_names_prev));
	
		plan_indicator = trues(size(plans_prev)[1]);
		for i in eachindex(plan_indicator)
			plan_indicator[i] = in(plans_prev[i,:pt_name],plan_year_names_prev);
		end
		plans_prev = plans_prev[findall(plan_indicator),:];
	end
	
	plan_indicator = trues(size(plans_pmt)[1]);
	plans_pmt[!,:pmt_name] = string.(plans_pmt[!,:plan_name],plans_pmt[!,:year],plans_pmt[!,:rating_area]);
	for i in eachindex(plan_indicator)
		plan_indicator[i] = in(plans_pmt[i,:pmt_name],plan_pmt_names);
	end
	plans_pmt = plans_pmt[findall(plan_indicator),:];
	sort!(plans_pmt);
	
	insurer_vector = convert(Array,plans[!,:Insurer]);
	insurers = unique(insurer_vector);
	
	household_pmt_indices = zeros(Int64,length(data_pmt_names));
	own_plan_year_matrix = zeros(Int64,length(plan_pmt_names),length(plan_pmt_names));
	plan_year_matrix = zeros(Int64,length(plan_pmt_names),length(plan_pmt_names));
	firm_plan_matrix = zeros(Int64,length(insurers),length(plan_pmt_names));
	for j in 1:length(plan_pmt_names)	
		firm = plans_pmt[j,:Insurer];
		plan_year = plans_pmt[j,:year];
		mkt = plans_pmt[j,:rating_area];
		own_plan_year_indices = intersect(findall(plans_pmt[!,:Insurer] .== firm),
			findall(plans_pmt[!,:year] .== plan_year));
		plan_year_indices = findall(plans_pmt[!,:year] .== plan_year);
		plan_indices = findall(data_pmt_names .== plan_pmt_names[j]);
			
		firm_plan_matrix[findall(insurers .== firm)[1],j] = 1;	
		own_plan_year_matrix[own_plan_year_indices,j] .= 1;
		plan_year_matrix[plan_year_indices,j] .= 1; 
		household_pmt_indices[plan_indices] .= j;
	end
	
	household_pt_indices = zeros(Int64,length(data_pt_names));
	for j in 1:length(plan_year_names)
		plan_indices = findall(data_pt_names .== plan_year_names[j]);
		household_pt_indices[plan_indices] .= j;
	end
	
	if year > 2014
		household_pt_indices_prev = zeros(Int64,length(data_pt_names_prev));
		for j in 1:length(plan_year_names_prev)
			plan_indices = findall(data_pt_names_prev .== plan_year_names_prev[j]);
			household_pt_indices_prev[plan_indices] .= j;
		end
	end
	
	year_plan_indices = zeros(Int64,size(plans_pmt)[1]);
	for j in 1:size(plans)[1]
		year_plan_indices[intersect(findall(plans_pmt[!,:plan_name] .== plans[j,:plan_name]),
			findall(plans_pmt[!,:year] .== plans[j,:year]))] .= j;
	end
	
	year_plan_names = string.(plans_pmt[!,:plan_name],plans_pmt[!,:year]);
	
	plan_foc_moment_cond_names = string.(plans_pmt[!,:plan_name],plans_pmt[!,:year]);
	foc_moment_cond_names = sort(unique(plan_foc_moment_cond_names));
	foc_moment_indicator = zeros(Int64,length(plan_pmt_names),length(foc_moment_cond_names));
	pmt_pt_mapping = zeros(Int64,length(plan_pmt_names));
	for m in 1:length(foc_moment_cond_names)
		plan_indices = findall(plan_foc_moment_cond_names .== foc_moment_cond_names[m]);
		foc_moment_indicator[plan_indices,m] .= 1;
		pmt_pt_mapping[plan_indices] .= m;
	end
	
	nonsilver_indicators = sort(cat(dims=1,findall(data[uninsured_plans .== 0,:AV] .== .6),findall(data[uninsured_plans .== 0,:AV] .== .8),findall(data[uninsured_plans .== 0,:AV] .== .9))); 
	if year > 2014
		silver_indicators_prev = sort(cat(dims=1,findall(data_prev[uninsured_plans_prev .== 0,:AV] .== .7),findall(data_prev[uninsured_plans_prev .== 0,:AV] .== .73),
			findall(data_prev[uninsured_plans_prev .== 0,:AV] .== .87),findall(data_prev[uninsured_plans_prev .== 0,:AV] .== .94))); 
	end
	silver_indicators = sort(cat(dims=1,findall(data[uninsured_plans .== 0,:AV] .== .7),findall(data[uninsured_plans .== 0,:AV] .== .73),
		findall(data[uninsured_plans .== 0,:AV] .== .87),findall(data[uninsured_plans .== 0,:AV] .== .94))); 	
	
	if random_parameters 
		
		if add_av_interactions
			random_variables = ["intercept","Anthem","Blue_Shield","Kaiser","Health_Net"];
		else
			random_variables = ["intercept","AV","Anthem","Blue_Shield","Kaiser","Health_Net"];
		end
		random_parameter_product_indices = cat(dims=1,findall(all_covariates .== "intercept")[1],findall(in(random_variables),product_covariates));
		random_product_covariates = string.(random_variables,"_random");
		
		random_parameter_indices = findall(in(random_product_covariates),all_covariates);
		if length(random_parameter_indices) == 1
			random_parameter_indices = random_parameter_indices[1];
		end
		
		# Assign a household reference number
		household_ref_numbers = zeros(maximum(households[!,:household_id]));
		household_ref_numbers[unique(households[!,:household_id])] = collect(1:1:length(unique(households[!,:household_id])));
		households[!,:ref_number] = Int64.(household_ref_numbers[households[!,:household_id]]);
		if year > 2014
			households_prev[!,:ref_number] = Int64.(household_ref_numbers[households_prev[!,:household_id]]);
		end
		
		#V = zeros(size(data)[1],length(random_variables));
		#for k in 1:length(random_variables)
		#	Random.seed!(random_parameter_product_indices[k]);
		#	Vi = randn(length(unique(households[!,:ref_number])));
		#	Vi = Vi[households[!,:ref_number]];
		#	V[:,k] = Vi[household_numbers_data];
		#end
		
		V = zeros(size(data)[1],length(random_variables),ns);
		if year > 2014
			V_prev = zeros(size(data_prev)[1],length(random_variables),ns);
		end
		
		for k in 1:length(random_variables)
			Random.seed!(random_parameter_product_indices[k]);
			Vi = randn(length(unique(households[!,:ref_number])),ns);
			
			Vi = Vi[households[!,:ref_number],:];
			V[:,k,:] = Vi[household_numbers_data,:];
			
			if year > 2014
				Vi_prev = Vi[households_prev[!,:ref_number],:];
				V_prev[:,k,:] = Vi_prev[household_numbers_data_prev,:];
			end
		end
		
	end
	
	plan_ra_factors = ones(length(plan_pmt_names)) * 0.6;
	plan_ra_factors[plans_pmt[!,:AV] .== 0.7] .= 0.7 * 1.03;
	plan_ra_factors[plans_pmt[!,:AV] .== 0.8] .= 0.8 * 1.08;
	plan_ra_factors[plans_pmt[!,:AV] .== 0.9] .= 0.9 * 1.15;

	metal_matrix = cat(dims=2,(plans_pmt[!,:AV] .== 0.6) .* 1,(plans_pmt[!,:AV] .== 0.7) .* 1,(plans_pmt[!,:AV] .== 0.8) .* 1,(plans_pmt[!,:AV] .== 0.9) .* 1);
	
	nu_for_riskadj = 0.86;
	
##### Create Data Objects
	# hdata object: contains household data (things that don't change)
	# pmt object: contains plan-market-year variables
	
	# X object: demand covariates
	# X_CLMS object: claims covariates (only plans available on exchange/in data)
	# X_RS: risk score covariates (only plans available on exchange/in data)
	
	# Field names
		
		hdata_fields = ["penalty","weight","pricing_factor","rating_factor","age_factor","ra_factor","GCF","alpha",
						"real_prob","real_ex_prob","risk_score","members","choice","real_slc_plan","real_unsubs_premium","real_instrumented_prob"];
		hvar_fields = ["subsidized_premium","subsidy","slc_plan","prob","ex_prob","risk_score"];
		pmt_fields = ["risk_score","ra_factor","real_demand","plan_pricing_weight","average_age_factor",
						"total_premiums","reins_factor","variable_admin","fixed_admin","base_premium"];
		partial_names = ["partial_plan_demand_partial_price","partial_firmrevenue_partial_price",
			"partial_firmtransfer_partial_price","partial_firmadmin_partial_price"];
	
		X = form_X_matrix(data,households,plans,product_covariates);
		if year > 2014
			X_prev = form_X_matrix(data_prev,households_prev,plans_prev,product_covariates);
		end
		hh_share_matrix = update_dem_share_matrix(households);
		
		hdata = initialize_household_data(X,hdata_fields,data,plans,household_numbers,min_beta);
		rs_share_variable_values = create_risk_score_share_variable_matrix(X,households);
		pmt = initialize_pmt_data(pmt_fields,plans_pmt,hdata); 
		X_RS = form_X_RS_matrix(plans_pmt,hdata,household_numbers,risk_score_covariates);
		X_CLMS = form_X_CLMS_matrix(plans_pmt,hdata,claims_moments_covariates,rs_parameter_estimates[!,2]);
		
		dhm_demands = create_dhm_demands(pmt);
		fixed_costs = pmt[:,findall(pmt_fields .== "fixed_admin")[1]] .* pmt[:,findall(pmt_fields .== "real_demand")[1]];
		average_plan_fixed_costs = (transpose(foc_moment_indicator) * fixed_costs) ./
			(transpose(foc_moment_indicator) * pmt[:,findall(pmt_fields .== "real_demand")[1]]);

	# Model parameters
	
		# Demand
		if nested
			lambda = last(min_beta);
			min_beta = min_beta[1:(length(min_beta)-1)];
		else
			lambda = 1;
		end
		
		# Claims and risk score
		
		if create_standard_errors_flag
			Random.seed!(sedraw + totalsedraws);
			adjust_variables = zeros(length(claims_moments_covariates));
			adjust_variables[1:(length(claims_moments_covariates) - length(rating_area_variables))] .= 1;
			#theta_CLMS = claims_parameter_estimates[!,2] .+ randn(length(claims_moments_covariates)) .* claims_parameter_estimates[!,3] .* adjust_variables;
			theta_CLMS = claims_parameter_estimates[!,2];
		else
			theta_CLMS = claims_parameter_estimates[!,2];
		end
		
		if create_standard_errors_flag
			rs_full_parameter_estimates = CSV.read(string("rs_reg_output_intermediate_avint.csv"),DataFrame, header=true);
			Random.seed!(sedraw + 2 * totalsedraws);
			#theta_RS = rs_full_parameter_estimates[!,2] .+ randn(length(risk_score_covariates)) .* rs_full_parameter_estimates[!,3];
			theta_RS = rs_parameter_estimates[!,2];
		elseif separate_metal_flag
			theta_RS_bronze = rs_estimates_bronze[!,2];
			theta_RS_silver = rs_estimates_silver[!,2];
			theta_RS_gold = rs_estimates_gold[!,2];
			theta_RS_platinum = rs_estimates_platinum[!,2];
			theta_RS = cat(dims=2,theta_RS_bronze,theta_RS_silver,theta_RS_gold,theta_RS_platinum);
		else
			theta_RS = rs_parameter_estimates[!,2];
		end
		theta_RS_share_indices = findall(in(risk_score_share_variables),risk_score_covariates);
		theta_RS_share_variables = theta_RS[theta_RS_share_indices];
	
	# Calculate Voucher
	
		voucher = get_voucher(hdata,X);
	
	# Initialize FOC residuals
		
		if !nochurn_flag & !zero_resid_flag
			foc_residuals = update_foc_values(plans[!,:Premium],pmt,hdata,X_RS,X_CLMS,zeros(length(plan_year_names)),false,false,false,false,false,false,false);	
			if (!network_flag & !eq_robust_flag)
				writedlm(string(year,"_foc_residuals.csv"),foc_residuals)
			end
		else
			foc_residuals = zeros(length(plan_names));
		end
		
		#if zero_resid_flag
		#	foc_residuals = zeros(length(plan_names));
		#end